pro hisa_images

cloudnum=[1000,1001]
glon=[31.409764,31.409764]
glat=[0.31116336,0.31116336]
vlsr=[96.8005,96.8005]
rad_as=2.354*[60.,60.]

restore,'/data13/nordhaus/gbt_data/ammonia/analysis/l033/l033_nh3_flag_1_nodup.sav';

b=sort(glon)
glon=glon[b]
glat=glat[b]
vlsr=vlsr[b]
rad_as=replicate(2.354*60.,n_e(glon))
cloudnum=replicate(1000,n_e(glon))+indgen(n_e(glon))

num=n_e(glon)

moleclinefilename=strarr(num)
hifilename=strarr(num)
vgpspath='/data13/nordhaus/VGPS/'
sgpspath='/data13/nordhaus/SGPS/'
grspath='/data13/nordhaus/GRS/'
for i=0,num-1 do begin
    ;something here for l009 field.
    if glon[i] ge 19.0 and glon[i] le 20.5 then moleclinefilename[i]=grspath+'grs-20-cube.fits'
    if glon[i] gt 20.5 and glon[i] le 21.5 then moleclinefilename[i]=grspath+'grs-21-cube.fits'
    if glon[i] gt 21.5 and glon[i] le 23.0 then moleclinefilename[i]=grspath+'grs-22-cube.fits'
    if glon[i] ge 31.0 and glon[i] le 32.5 then moleclinefilename[i]=grspath+'grs-32-cube.fits'
    if glon[i] gt 32.5 and glon[i] le 33.5 then moleclinefilename[i]=grspath+'grs-33-cube.fits'
    if glon[i] gt 33.5 and glon[i] le 35.0 then moleclinefilename[i]=grspath+'grs-34-cube.fits'
    if glon[i] ge 51.0 and glon[i] le 52.5 then moleclinefilename[i]=grspath+'grs-52-cube.fits'
    if glon[i] gt 52.5 and glon[i] le 53.5 then moleclinefilename[i]=grspath+'grs-53-cube.fits'
    if glon[i] gt 53.5 and glon[i] le 54.5 then moleclinefilename[i]=grspath+'grs-54-cube.fits'
    if glon[i] gt 54.5 and glon[i] le 56.0 then moleclinefilename[i]=grspath+'grs-55-cube.fits'

    if glon[i] ge 7.3 and glon[i] le 10.2 then hifilename[i]=sgpspath+'G8.750.0.0.HI.sgps.fits'
    if glon[i] ge 18.5 and glon[i] le 23.5 then hifilename[i]=vgpspath+'MOS_021.Tb.fits'
    if glon[i] gt 30.5 and glon[i] le 35.5 then hifilename[i]=vgpspath+'MOS_033.Tb.fits'
    if glon[i] gt 50.5 and glon[i] le 55.0 then hifilename[i]=vgpspath+'MOS_053.Tb.fits'
    if glon[i] gt 55.0 and glon[i] le 59.5 then hifilename[i]=vgpspath+'MOS_057.Tb.fits'
end

;determine BGPS filename
bgpsfilename=strarr(num)
bgpspath='/data13/nordhaus/bolocam_data/releases/finalv1.0/map_original/'
;load in bgps_bounds file
readcol,'/data13/nordhaus/bolocam_data/releases/finalv1.0/bgps_bounds.txt',imagename,glon_min,glon_max,glat_min,glat_max,delimiter=' ',format='(A,D,D,D,D)',/silent
for j=0,n_e(imagename)-1 do begin
    w=where(glon ge glon_min[j] and glon lt glon_max[j] and glat ge glat_min[j] and glat lt glat_max[j],ct)
    if ct gt 0 then bgpsfilename[w]=bgpspath+imagename[j]
endfor



crap=strarr(1)

oldhifilename='junk'
oldmoleclinefilename='junk'
oldbgpsfilename='junk'
;loop over all sources
for i=0,10 do begin
    ;define boxsize and cornor coordinates for image
    halfboxsize=2.*rad_as[i]/3600.
    x1=glon[i]+halfboxsize
    x2=glon[i]-halfboxsize
    y1=glat[i]-halfboxsize
    y2=glat[i]+halfboxsize
    ;--------------------------------------------------
    ;HI map set up
    ;check if we've already read in the fits file
    if hifilename[i] ne oldhifilename then begin
        ;read in HI
        largehi=readfits(hifilename[i],largehihdr)
        extast,largehihdr,largehiastr
    endif

    ;must pull out velocity array before hextract since hextract requires 2d array
    largehiv = (lindgen(sxpar(largehihdr,'NAXIS3'))-(sxpar(largehihdr,'CRPIX3')-1.))*sxpar(largehihdr,'CDELT3')+sxpar(largehihdr,'CRVAL3')
    largehiw=where(largehiv-sxpar(largehihdr,'CDELT3')/2. gt vlsr[i]*1.d3 and largehiv+sxpar(largehihdr,'CDELT3')/2. le vlsr[i]*1.d3)
    temphi=largehi[*,*,largehiw[0]]
    sxaddpar,largehihdr,'NAXIS',2 ;hextract can only take 2d
    sxdelpar,largehihdr,'NAXIS3'  ;hextract can only take 2d
    sxdelpar,largehihdr,'NAXIS4'  ;hextract can only take 2d
    
    ;get pixel numbers for box corners for hi image
    ad2xy,x1,y1,largehiastr,hix1,hiy1
    ad2xy,x2,y2,largehiastr,hix2,hiy2
    ;do extraction on tempfile that has only l and b dimensions
    hextract,temphi,largehihdr,hi,hihdr,hix1,hix2,hiy1,hiy2
            ;writefits,'HI_'+strtrim(string(cloudnum[i]),2)+'.fits',hi,hihdr
    extast,hihdr,hiastrhdr
    
    hicrpix1 = sxpar(hihdr,'CRPIX1')-1. ;-1 to correct for difference between numbering in header vs idl
    hicrpix2 = sxpar(hihdr,'CRPIX2')-1.
    hicrval1 = sxpar(hihdr,'CRVAL1')
    hicrval2 = sxpar(hihdr,'CRVAL2')
    hicdelt1 = sxpar(hihdr,'CDELT1')
    hicdelt2 = sxpar(hihdr,'CDELT2')

    hix = lindgen(sxpar(hihdr,'NAXIS1'))
    hiy = lindgen(sxpar(hihdr,'NAXIS2'))
    hil = (hix-hicrpix1)*hicdelt1+hicrval1
    hib = (hiy-hicrpix2)*hicdelt2+hicrval2

    ;--------------------------------------------------
    ;MOLECLINE map set up
    ;check if we've already read in the fits file
    if moleclinefilename[i] ne oldmoleclinefilename then begin
        ;read in MOLECLINE
        largemolecline=readfits(moleclinefilename[i],largemoleclinehdr)
        extast,largemoleclinehdr,largemoleclineastr
    endif

    ;must pull out velocity array before hextract since hextract requires 2d array
    largemoleclinev = (lindgen(sxpar(largemoleclinehdr,'NAXIS3'))-(sxpar(largemoleclinehdr,'CRPIX3')-1.))*sxpar(largemoleclinehdr,'CDELT3')+sxpar(largemoleclinehdr,'CRVAL3')
    largemoleclinew=where((largemoleclinev-sxpar(largemoleclinehdr,'CDELT3')/2.) lt (vlsr[i]*1.d3) and (largemoleclinev+sxpar(largemoleclinehdr,'CDELT3')/2.) ge (vlsr[i]*1.d3))
    tempmolecline=largemolecline[*,*,largemoleclinew[0]]
    sxaddpar,largemoleclinehdr,'NAXIS',2 ;hextract can only take 2d
    sxdelpar,largemoleclinehdr,'NAXIS3'  ;hextract can only take 2d
    sxdelpar,largemoleclinehdr,'NAXIS4'  ;hextract can only take 2d
    
    ;get pixel numbers for box corners for molecline image
    ad2xy,x1,y1,largemoleclineastr,moleclinex1,molecliney1
    ad2xy,x2,y2,largemoleclineastr,moleclinex2,molecliney2
    ;do extraction on tempfile that has only l and b dimensions
    hextract,tempmolecline,largemoleclinehdr,molecline,moleclinehdr,moleclinex1,moleclinex2,molecliney1,molecliney2
            ;writefits,'MOLECLINE_'+strtrim(string(cloudnum[i]),2)+'.fits',molecline,moleclinehdr
    extast,moleclinehdr,moleclineastrhdr
    
    moleclinecrpix1 = sxpar(moleclinehdr,'CRPIX1')-1. ;-1 to correct for difference between numbering in header vs idl
    moleclinecrpix2 = sxpar(moleclinehdr,'CRPIX2')-1.
    moleclinecrval1 = sxpar(moleclinehdr,'CRVAL1')
    moleclinecrval2 = sxpar(moleclinehdr,'CRVAL2')
    moleclinecdelt1 = sxpar(moleclinehdr,'CDELT1')
    moleclinecdelt2 = sxpar(moleclinehdr,'CDELT2')

    moleclinex = lindgen(sxpar(moleclinehdr,'NAXIS1'))
    molecliney = lindgen(sxpar(moleclinehdr,'NAXIS2'))
    moleclinel = (moleclinex-moleclinecrpix1)*moleclinecdelt1+moleclinecrval1
    moleclineb = (molecliney-moleclinecrpix2)*moleclinecdelt2+moleclinecrval2

    ;get molecline image ready to plot as contour on HI image
    moleclinecontour_nx=sxpar(moleclinehdr,'NAXIS1')
    moleclinecontour_ny=sxpar(moleclinehdr,'NAXIS2')
    moleclinecontour_x = dindgen(moleclinecontour_nx*moleclinecontour_ny)
    moleclinecontour_y = dindgen(moleclinecontour_nx*moleclinecontour_ny)
    moleclinecontour_a = dindgen(moleclinecontour_nx*moleclinecontour_ny)
    moleclinecontour_d = dindgen(moleclinecontour_nx*moleclinecontour_ny)
    moleclinecontour_image1d = dblarr(moleclinecontour_nx*moleclinecontour_ny)
    step = long(0)
    for k=0,moleclinecontour_nx-1 do begin
        for j=0,moleclinecontour_ny-1 do begin
            moleclinecontour_x[step]=k
            moleclinecontour_y[step]=j
            xy2ad,moleclinecontour_x[step],moleclinecontour_y[step],moleclineastrhdr,temp1,temp2
            moleclinecontour_a[step]=temp1
            moleclinecontour_d[step]=temp2
            moleclinecontour_image1d[step]=molecline[k,j]
            step=step+1
        endfor
    endfor
    moleclinecontour_sigma=max(molecline)/10.
    moleclinecontour_levels=moleclinecontour_sigma*[4.,6.,9.]
    moleclinecontour_thick=2.0*[1.0,1.0,1.0]
    moleclinecontour_line=0*[1,1,1]
    
    ;--------------------------------------------------
    ;BGPS map set up
    ;check if we've already read in the fits file
    if bgpsfilename[i] ne oldbgpsfilename then begin
        ;read in BGPS
        largebgps=readfits(bgpsfilename[i],largebgpshdr)
        extast,largebgpshdr,largebgpsastr
    endif

    ;;must pull out velocity array before hextract since hextract requires 2d array
    ;largebgpsv = (lindgen(sxpar(largebgpshdr,'NAXIS3'))-(sxpar(largebgpshdr,'CRPIX3')-1.))*sxpar(largebgpshdr,'CDELT3')+sxpar(largebgpshdr,'CRVAL3')
    ;largebgpsw=where(largebgpsv-sxpar(largebgpshdr,'CDELT3')/2. gt vlsr[i]*1.d3 and largebgpsv+sxpar(largebgpshdr,'CDELT3')/2. le vlsr[i]*1.d3)
    ;tempbgps=largebgps[*,*,largebgpsw[0]]
    ;sxaddpar,largebgpshdr,'NAXIS',2 ;hextract can only take 2d
    ;sxdelpar,largebgpshdr,'NAXIS3'  ;hextract can only take 2d
    ;sxdelpar,largebgpshdr,'NAXIS4'  ;hextract can only take 2d
    
    ;get pixel numbers for box corners for bgps image
    ad2xy,x1,y1,largebgpsastr,bgpsx1,bgpsy1
    ad2xy,x2,y2,largebgpsastr,bgpsx2,bgpsy2
    ;do extraction on tempfile that has only l and b dimensions
    hextract,largebgps,largebgpshdr,bgps,bgpshdr,bgpsx1,bgpsx2,bgpsy1,bgpsy2
            ;writefits,'BGPS_'+strtrim(string(cloudnum[i]),2)+'.fits',bgps,bgpshdr
    extast,bgpshdr,bgpsastrhdr
    
    bgpscrpix1 = sxpar(bgpshdr,'CRPIX1')-1. ;-1 to correct for difference between numbering in header vs idl
    bgpscrpix2 = sxpar(bgpshdr,'CRPIX2')-1.
    bgpscrval1 = sxpar(bgpshdr,'CRVAL1')
    bgpscrval2 = sxpar(bgpshdr,'CRVAL2')
    bgpscdelt1 = sxpar(bgpshdr,'CD1_1')
    bgpscdelt2 = sxpar(bgpshdr,'CD2_2')

    bgpsx = lindgen(sxpar(bgpshdr,'NAXIS1'))
    bgpsy = lindgen(sxpar(bgpshdr,'NAXIS2'))
    bgpsl = (bgpsx-bgpscrpix1)*bgpscdelt1+bgpscrval1
    bgpsb = (bgpsy-bgpscrpix2)*bgpscdelt2+bgpscrval2

    ;get bgps image ready to plot as contour on HI image
    bgpscontour_nx=sxpar(bgpshdr,'NAXIS1')
    bgpscontour_ny=sxpar(bgpshdr,'NAXIS2')
    bgpscontour_x = dindgen(bgpscontour_nx*bgpscontour_ny)
    bgpscontour_y = dindgen(bgpscontour_nx*bgpscontour_ny)
    bgpscontour_a = dindgen(bgpscontour_nx*bgpscontour_ny)
    bgpscontour_d = dindgen(bgpscontour_nx*bgpscontour_ny)
    bgpscontour_image1d = dblarr(bgpscontour_nx*bgpscontour_ny)
    step = long(0)
    for k=0,bgpscontour_nx-1 do begin
        for j=0,bgpscontour_ny-1 do begin
            bgpscontour_x[step]=k
            bgpscontour_y[step]=j
            xy2ad,bgpscontour_x[step],bgpscontour_y[step],bgpsastrhdr,temp1,temp2
            bgpscontour_a[step]=temp1
            bgpscontour_d[step]=temp2
            bgpscontour_image1d[step]=bgps[k,j]
            step=step+1
        endfor
    endfor
    bgpscontour_sigma=max(bgps)/10.
    if bgpscontour_sigma ge 0.2 then bgpscontour_levels=bgpscontour_sigma*[0.5,0.75,1.]
    if bgpscontour_sigma lt 0.2 then bgpscontour_levels=bgpscontour_sigma*[0.8,0.9,1.]
    bgpscontour_thick=2.0*[1.0,1.0,1.0]
    bgpscontour_line=5*[1,1,1]
    




pltfile='images/HI_map_'+strtrim(string(cloudnum[i]+i),2)+'.ps'
!p.multi=[0,1,1,0,0]
!x.thick=3.0
!y.thick=3.0
!p.thick=3.0
!p.charsize=1.25
!p.charthick=3
set_plot,'ps'
device,/portrait,yoffset=4.5,xoffset=1.5,filename=pltfile,xsize=22,ysize=20
cgLoadct,0

imdisp,hi,/axis,xrange=[max(hil),min(hil)],yrange=[min(hib),max(hib)],range=[min(hi),max(hi)],xtitle='Galactic Longitude (deg)',ytitle='Galactic Latitude (deg)',title=strtrim(string(cloudnum[i]),2),ticklen=.02,true=1,erase=0,_extra=_extra,margin=0.075;,negative=1
contour,moleclinecontour_image1d,moleclinecontour_a,moleclinecontour_d,/irregular,/overplot,levels=moleclinecontour_levels,c_thick=moleclinecontour_thick,c_line=moleclinecontour_line
contour,bgpscontour_image1d,bgpscontour_a,bgpscontour_d,/irregular,/overplot,levels=bgpscontour_levels,c_thick=bgpscontour_thick,c_line=bgpscontour_line

device,/close_file
spawn,'gv images/HI_map_'+strtrim(string(cloudnum[i]+i),2)+'.ps &'


endfor


end




